---
title: "R1.Figures"
output: pdf_document
---

# Load Packages 
```{r}
library(tidyverse)
library(lubridate)
library(sf)
library(tmap)
```

# Load Data
```{r}
load("Data/data_final.Rdata")
load("Data/zira_portugal.Rdata")
load("Data/gnr_panel.Rdata")
load("Data/gnr_events.Rdata") 
load("Data/ideo_list.Rdata")
load("Data/DG_Z_Shp.Rdata")
```

# Generate City Coordinates
```{r}
cities <-data.frame(City = c("Lisbon", "Évora", "Beja", "Portalegre", "Setúbal"), 
coords = c("38.722158, -9.139347",
           "38.570620, -7.913413",
           "38.015333, -7.863233",
           "39.291186, -7.430803",
           "38.524769, -8.892533")) %>% 
  separate(coords, into = c("lat", "lon"), sep = ", ") %>% 
  st_as_sf(coords = c("lon", "lat"), agr = "constant", crs = 4326, stringsAsFactors = FALSE, remove = TRUE)
```

# Map of Points Index (Not Included in Paper or Appendix)
```{r}
# map_points <- tm_shape(DG_Z_Shp %>% left_join(data_final)) +
#   tm_polygons("CA1_CAR1U_MeanF", title = "Points Index", palette = "Greys", colorNA = "white", legend.hist = T) + tm_legend(legend.outside=T) +
#    tm_shape(cities) + tm_symbols(shape = "City", col = "black",  border.col = "white", size = 0.5, alpha = 0.7)
# tmap_save(map_points, "Figures/Map_PointsIndex.pdf")
```

# Map of ZIRA (Figure 1a)
```{r}
map_zira_pt <- tm_shape(zira_portugal[[2]] %>% st_transform(crs = st_crs(zira_portugal[[1]]))) + tm_borders() +
  tm_shape(zira_portugal[[1]]) + tm_polygons(fill = "grey")

tmap_save(map_zira_pt, "Figures/fg1a_Map_ZIRA.pdf")
```

# Map of Potential Productivity (Figure 1b)
```{r}
map_baseline <- tm_shape(DG_Z_Shp %>% left_join(data_final)) +
  tm_polygons("CaloricMean", title = "Potential Productivity \n (Caloric Output)", palette = "Greys", colorNA = "white", legend.hist = T) + tm_legend(legend.outside=T) + 
   tm_shape(cities) + tm_symbols(shape = "City", col = "black",  border.col = "white", size = 0.5, alpha = 0.7)
tmap_save(map_baseline, "Figures/fg1b_Map_ProductPotent.pdf")
```

# Map of Mismatch Index (Figure 1c)
```{r}
map_resid <- tm_shape(DG_Z_Shp %>% left_join(data_final)) +
  tm_polygons("resid", title = "Point Discrepancy", style = "fixed", breaks = c(-30, -20, -10, 0, 10, 20, 30), 
              legend.hist = T, palette = "RdGy", colorNA = "white") + tm_legend(legend.outside=T) + 
 tm_layout(legend.hist.width = 5) +
  tm_shape(cities) + tm_symbols(shape = "City", col = "black",  border.col = "white", size = 0.5, alpha = 0.7)
tmap_save(map_resid, "Figures/fg1c_Map_Mismatch.pdf")
```

# Map of Expropriations (Figure 2)
```{r}
map_exprop <- tm_shape(DG_Z_Shp %>% left_join(data_final) %>% filter(!is.na(pr_ha_ord))) +
  tm_polygons("pr_ha_ord", style = "fixed", 
              breaks = c(0, 0.001, 10, 25, 50, 75, 100), 
              labels = c("0%", "0-10%", "10-25%", "25-50%", "50-75%", "75-100%"), 
              palette = c("white", "grey90", "grey70", "grey50", "grey30", "grey10"), lines = "black",
              title = "% Hectares Expropriated", colorNA = "red") +  tm_legend(legend.outside=T) +
  tm_shape(cities) + tm_symbols(shape = "City", col = "black",  border.col = "white", size = 0.5, alpha = 0.7)
tmap_save(map_exprop, "Figures/fg2_Map_Expropriation.pdf")
```

# Map of Treatment Indicator (Figure 3b)
```{r}
map_tx75 <- tm_shape(DG_Z_Shp %>% left_join(data_final) %>% data.frame() %>% st_as_sf()) +
  tm_polygons("tx_Hi", style = "cat", palette = c("white", "black"), colorNA = "grey", title = "Treatment (75th Percentile)") + tm_legend(legend.outside=T)
tmap_save(map_tx75, "Figures/fg3b_Map_Treatment.pdf")
```

# Timeline (Figure 5)
```{r}

raw <- data.frame(date = c(as.Date("1973-01-01"),
                              as.Date("1974-04-25"),
                              as.Date("1975-04-25"),
                              as.Date("1976-04-25"),
                              as.Date("1979-12-02"),
                              as.Date("1980-10-05"),
                              as.Date("1983-04-25"),
                              as.Date("1985-10-06"),
                              as.Date("1987-07-19"),
                              as.Date("1991-10-06"),
                              as.Date("1995-10-01"),
                              as.Date("1999-10-10"),
                              as.Date("2002-03-17"),
                              as.Date("2005-02-20"),
                              as.Date("2009-09-27"),
                              as.Date("2011-06-05"), 
                              as.Date("1975-07-29"), 
                              as.Date("1977-09-29")), 

event = c("Estado Novo",
          "Carnation Revolution",
          "Const. Assembly Election",
          "1976 Legislative Election",
          "1979 Legislative Election",
          "1980 Legislative Election",
          "1983 Legislative Election",
          "1985 Legislative Election",
          "1987 Legislative Election",
          "1991 Legislative Election",
          "1995 Legislative Election",
          "1999 Legislative Election",
          "2002 Legislative Election",
          "2005 Legislative Election",
          "2009 Legislative Election",
          "2011 Legislative Election", 
          "Land Reform Law",
          "First Counter-Reform Law"), 

period = c("Estado Novo",
           "MFA",
          "Transition",
          "PS",
          "PSD",
          "PSD",
          "PS + PSD",
          "PSD",
          "PSD",
          "PSD",
          "PS",
          "PS",
          "PSD + CDS",
          "PS",
          "PS", 
          "PSD", 
          NA, 
          NA))

lr_time <- data.frame(event = "Land Reform", start = as.Date("1975-07-01"), end = as.Date("1976-08-01"))

periods <- raw %>% slice(1:(nrow(.)-4)) %>% rename(start = date) %>% mutate(end = lead(start)) %>% 
  filter(start < as.Date("1980-10-05"))

events <- raw %>% filter(date > as.Date("1973-01-01")) %>% 
  filter(date < as.Date("1980-10-05")) %>% arrange(date)

ggplot() +
  geom_vline(data = events %>% filter(date > as.Date("1970-01-01")), aes(xintercept = date), col = "red", alpha = 0.5, size = 0.5) +
  geom_rect(data = periods, aes(xmin = start, xmax = end, ymin = -0.5, ymax = 0.5), fill= "grey", color="black", show.legend = FALSE) +
  geom_text(data = periods, aes(x = start + floor((end - start)/2), label = period, y = 0), size = 3) +
  geom_text(data = events[c(1,4,6),], aes(x = date, label = event, y = -0.7), nudge_y = 0, size = 2) +
  geom_text(data = events[c(2,5),], aes(x = date, label = event, y = -0.7), nudge_y = -0.2, size = 2) +
  geom_text(data = events[c(3),], aes(x = date, label = event, y = -0.7), nudge_y = -0.4, size = 2) +
  geom_segment(data = lr_time, aes(x = start, xend = end, y = 0.6, yend = 0.6, lineend = "butt"), col = "black") +
  geom_text(data = lr_time, aes(x = start + floor((end - start)/2), label = event, y = 0.75), size = 3) +
  theme(legend.position="none") + 
  ylim(-1.2, 1) +
  xlim(as.Date("1973-01-01"), as.Date("1980-10-05")) +
  theme_light() +
  xlab("Year") + ylab("Hectares of Land Expropriated and Retained") +
  theme(axis.text.y = element_text(angle = 45)) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())  + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank())

ggsave(plot = last_plot(), "Figures/fg5_Timeline.pdf", width = 6, height = 2)

```
# GNR Events Plot
```{r}
gnr_events %>% mutate(date = ymd(date), mo_yr = round_date(date, "month")) %>% group_by(mo_yr) %>% summarise(events = sum(gnr_rolling_Adm23_Z)) %>%
  ggplot() + geom_col(aes(x=mo_yr, y=events)) + theme_light() + ylab("GNR Interventions") + xlab("Time in Months")

ggsave(plot = last_plot(), "Figures/Afg1a1.1_GNREvents.pdf")
```

# Map of GNR Events over Time (App. Figure 1ai)
```{r}
gnr_period_map <- gnr_panel %>% filter(pre_election != "AC75") %>% left_join(DG_Z_Shp %>% st_make_valid()) %>% st_as_sf %>%
  tm_shape() + tm_polygons("gnr_periods_23_Z", palette = "Greys", title = "GNR Events Prior \n To Each Election", 
              breaks = c(0, 0.001, 5, 10, 15, 20, 25, 30, 35, 40), labels = c("0", "1-5", "5-10", "10-15", "15-20", "20-25", "25-30", "30-35", "35-40"), ) + 
  tm_facets("pre_election", free.coords = F)
tmap_save(gnr_period_map, "Figures/Afg1a1.2_Map_GNRWaves.pdf")
```

# Map of PCP Office Indicator (App. Figure 1aii)
```{r}
map_PCPOffice <- tm_shape(DG_Z_Shp %>% st_make_valid() %>% left_join(data_final) %>% data.frame() %>% st_as_sf() %>% 
                            mutate(PCPOffices_Map = ifelse(PCPOffices_Map >= 1, 1, 0))) +
  tm_polygons("PCPOffices_Map", style = "cat", palette = c("white", "black"), colorNA = "grey", title = "Pre-LR PCP Office") +
  tm_shape(cities) + tm_symbols(shape = "City", col = "black",  border.col = "white", size = 0.5, alpha = 0.7) + 
  tm_legend(legend.outside=T)
tmap_save(map_PCPOffice, "Figures/Afg1a2_Map_PCPOffice.pdf")
```

# Map of Vote Share Change 1975-76, 1975-79 (App. Figure 1dii)
```{r}
early_changes <- do.call(rbind, ideo_list[2:3]) %>% rename(new = votes) %>% 
  left_join(ideo_list[[1]] %>% rename(base = votes) %>% select(-election)) %>% 
  mutate(first_diff = new - base) %>% 
  left_join(data_final) %>% 
  rename(Party = ideology2) %>% 
  mutate(election = as.character(election)) %>% 
  filter(!is.na(exp_cat2)) %>% 
  mutate(Party = ifelse(Party == "Center", "Center-Left", Party), 
         Party = ifelse(Party == "Left", "Far Left", Party)) %>% 
  left_join(DG_Z_Shp %>% st_make_valid()) %>% st_as_sf %>% tm_shape() + tm_fill("first_diff", title = "Vote Change", n=10) + tm_facets(c("election", "Party"))

tmap_save(early_changes, "Figures/Afg1d2.1_Map_EarlyChanges.pdf", width = 10, height = 7)
```

# Map of Casa do Povo Creation Year (App. Figure 3ciii)
```{r}
# Missing: 7. Evora, 14. Santarem
missing_D <- DG_Z_Shp %>% 
  filter((str_detect(DICOFRE_Z, "^7") | str_detect(DICOFRE_Z, "^14"))) %>%
  separate(DICOFRE_Z, into = c("D", "MF"), -4, remove = F) %>% 
  group_by(D) %>% summarise() %>%
  ungroup() %>% st_as_sf() %>% st_buffer(0.001)
cdp_map<-tm_shape(DG_Z_Shp %>% st_make_valid() %>% left_join(data_final)) + tm_polygons("cdp_year", n = 10, 
                                                      palette = "Greys",
                                                      colorNA = "white",
              title = "Year of Casa do Povo creation") +  tm_legend(legend.outside=T)  +
  tm_shape(missing_D) + tm_borders(col = "red")
tmap_save(cdp_map, "Figures/Afg3c3_Map_CdPYear.pdf")
```

# Map of Land Reform Intensity (App. Figure 4c)
```{r}
lr_intensity <- DG_Z_Shp %>% st_make_valid %>% left_join(data_final) %>% st_as_sf %>% 
  tm_shape() + tm_polygons("exp_cat2", style = "cat", 
                           palette = c("grey90", "grey50", "grey10"),
                           title = "Land Reform\nIntensity/Proximity",
                           labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"), colorNA = "white",
                           ) +
  tm_legend(legend.outside=T) 
tmap_save(lr_intensity, "Figures/Afg4c_Map_LRIntensity.pdf")
```

# Map of Land Returns (App. Figure 8bii)
```{r}
all_returned <- DG_Z_Shp %>% st_make_valid() %>% 
  left_join(data_final %>% mutate(all_returned = ifelse(!(is.na(pr_haret_exp_coop_85) | pr_haret_exp_coop_85 < 100), 1, 0))) %>% st_as_sf() %>% 
  tm_shape() + tm_polygons("all_returned", style = "cat", palette = c("grey90", "grey10"),, title = "All Land Returned", colorNA = "white", lines = "black") +
  tm_legend(legend.outside=T) 
tmap_save(all_returned, "Figures/Afg8b2_Map_AllReturned.pdf")
```
